Resonance Raman enhancement by the intralayer and interlayer electron–phonon processes in twisted bilayer graphene

Twisted bilayer graphene is a fascinating system due to the possibility of tuning the electronic and optical properties by controlling the twisting angle \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta$$\end{document}θ between the layers. The coupling between the Dirac cones of the two graphene layers gives rise to van Hove singularities (vHs) in the density of electronic states, whose energies vary with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta$$\end{document}θ. Raman spectroscopy is a fundamental tool to study twisted bilayer graphene (TBG) systems since the Raman response is hugely enhanced when the photons are in resonance with transition between vHs and new peaks appear in the Raman spectra due to phonons within the interior of the Brillouin zone of graphene that are activated by the Moiré superlattice. It was recently shown that these new peaks can be activated by the intralayer and the interlayer electron–phonon processes. In this work we study how each one of these processes enhances the intensities of the peaks coming from the acoustic and optical phonon branches of graphene. Resonance Raman measurements, performed in many different TBG samples with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta$$\end{document}θ between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4^{\circ }$$\end{document}4∘ and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$16^{\circ }$$\end{document}16∘ and using several different laser excitation energies in the near-infrared (NIR) and visible ranges (1.39–2.71 eV), reveal the distinct enhancement of the different phonons of graphene by the intralayer and interlayer processes. Experimental results are nicely explained by theoretical calculations of the double-resonance Raman intensity in graphene by imposing the momentum conservation rules for the intralayer and the interlayer electron–phonon resonant conditions in TBGs. Our results show that the resonant enhancement of the Raman response in all cases is affected by the quantum interference effect and the symmetry requirements of the double resonance Raman process in graphene.

www.nature.com/scientificreports/ spectra associated with phonons in the interior of Brillouin zone (BZ) of monolayer graphene that are activated by the Moiré pattern in the spectra of TBGs [20][21][22][23][24][25][26] . Other Raman bands are also observed in graphene systems due to the double-resonance (DR) mechanism, such as the 2D band and the disorder-induced D band, being the dependence of the band positions and intensities on the laser excitation energy their main characteristics 17,27 . Since the first report by Gupta et al. in 2010 20 several works in the literature have also reported the appearance of new peaks in the Raman spectra of TBGs, and have shown that their frequencies depend on the twisting angle θ [21][22][23][24][25][26] . These peaks come from phonons within the interior of the Brillouin zone (BZ) of monolayer graphene that have the wavevector q M of the Moiré superlattice reciprocal space. They are folded to the centre of the reduced BZ of the Moiré superlattice and become Raman active. It was shown recently that there are two different mechanisms that enhance the Moiré peaks in the TBG Raman spectra 28 . In one case, the enhancement of these peaks occurs in the same range of excitation energies where the G band is hugely enhanced, and, in the other case, the enhancement occurs for excitation laser energies where the G band is not enhanced. These results were explained considering the existence of two resonance mechanisms: the intralayer and interlayer electron-phonon (el-ph) interaction processes. The intralayer process involves electronic scattering between states in the Dirac cone of the same layer and in the interlayer case the scattering occurs between states from different layers 28 .
In the present work, we study experimentally and theoretically the enhancement of the Moiré peaks of TBGs coming from the acoustic and optical phonon branches of graphene by both the intralayer and the interlayer el-ph processes. We present experimental results obtained from many different TBG samples with small and intermediate twisting angles (between 4.3 • and 16.3 • ) and using several different laser excitation energies in the visible and NIR ranges (between 1.39 eV and 2.71 eV). Our results show that the enhancement of Moiré peaks from the different acoustic and optical branches of graphene by the intralayer and the interlayer electron-phonon processes are distinct. The experimental results are explained by calculations based on double-resonance Raman processes in graphene 7,27 by imposing the momentum conservation for the intralayer and interlayer el-ph processes. We also show that the positions of the Moiré peaks from the acoustic phonon branches can be used to measure accurately the twisting angle θ of a TBG sample.

Results and discussion
Acoustic and optical phonons of graphene in TBG Raman spectra. Figure 1a shows the Raman spectra of 10 different samples of TBG with twisting angles between 4.3 • and 10.2 • , which have been measured with different laser excitation energies in the NIR range: 1.39 eV (890 nm), 1.50 eV (827 nm), 1.58 eV (785 nm), 1.70 eV (728 nm) and 1.82 eV (680 nm). The spectra depicted in Fig. 1a correspond to the excitation energy that gives rise to the strongest resonance Raman effect for each sample. We can observe that the spectra can be separated in three spectral ranges: (i) between 100 cm −1 and 600 cm −1 we observe the peaks associated with the acoustic phonon branches (ZA, TA and LA), (ii) around 850 cm −1 we observe the peaks associated with the out-of-plane transverse optical (oTO) phonon mode and (iii) in the range 1400-1700 cm −1 we observe the G Figure 1. (a) Raman spectra of ten TBG samples with different angles between 4.3 • and 10.2 • . The spectra were acquired using four different laser lines, represented by different colours in the figure, but only those presenting the strongest resonance effect are displayed in the figure. The spectra is shown in three regions: (i) the acoustic modes (ZA, TA and LA), (ii) the out-of-plane TO phonon (oTO) and (iii) the optical modes close to the G band (iTO and LO). All spectra with θ ≥ 7.4 • were multiplied by 10 × in the acoustic and the oTO regions for a better visualisation. (b) Positions of all peaks extracted from all Raman spectra plotted as a function of θ . The twisting angles were determined by the LA and TA phonon branches according to dispersion curves (full black lines) obtained by Cong et al. 29 (see text). The dashed black lines are theoretical curves of the TBG phonon frequencies 27 as a function of the twisting angle 28 . The blue (green) stripe highlights the range of angles of the TBGs investigated in this work that are activated by the intralayer (interlayer) el-ph process. www.nature.com/scientificreports/ band (1580 cm −1 ) and the peaks originated from the in-plane transverse optical (iTO) phonon branch and the longitudinal optical (LO) phonon branch, respectively below and above the G band. The twisting angles θ were initially estimated from the optical images of the TBG samples, where we can observe the crystalline edges of the two graphene layers, as previously reported by Ribeiro et al. 19 . However, the angle determined using this method is not accurate and the uncertainty is of the order of 1 • . The accuracy in the angles determination was improved by considering the dependence of acoustic phonon frequencies (TA and LA) on the twisting angle θ , since the dispersive behavior of these branches provides the determination of θ with an accuracy of the order of 0.1 • .
Cong et al. 29 have probed the LA and TA phonon branches of monolayer graphene near the Ŵ point by double resonant Raman scattering of combination phonon modes and the overtone 2 D ′ mode. They obtained the values of 12.9 km s −1 and 19.9 km s −1 for the sound velocities of TA and LA branches, respectively. For TBG, the twisting angle may be obtained by the expression sin(θ/2) = f √ 3a/4v , being a the lattice parameter of monolayer graphene, f the frequency and v the sound velocity of the TA or LA branch. The twisting angle θ is then obtained from the experimental frequencies of the TA and LA phonons using the solid lines in Fig. 1b, where the TA peak position is related to the twisting angle by the expression f TA = 4040 sin(θ/2) cm −1 while the LA curve is given by f LA = 6232 sin(θ/2) cm −1 . In the case of a TBG sample with the magic angle of 1.1 • , the TA and LA peaks are expected to appear in the Raman spectrum around 39 cm −1 and 60 cm −1 , respectively. However, as we will show below, the enhancement of these peaks is only expected to be observed for excitation energies in the mid-IR range, below the existent laser energies in the NIR range. Figure 1b shows the positions of all Raman peaks of the investigated TBG samples as a function of the twisting angle θ . The dispersion relations for the ZA, oTO, iTO and LO phonons are also plotted in Fig. 1b as dashed curves and were calculated by means of density functional theory as discussed in Ref. 27 . The precise determination of the ZA peak positions is affected by the fact that the spectrometer cuts the scattered light below 100 cm −1 . Therefore, we only plotted the data of the highest angle samples, where the peak position is above 100 cm −1 .

Resonance Raman enhancement by the intralayer and interlayer el-ph processes.
We can clearly distinguish in Fig. 1a two different resonance behaviors by analysing the intensity of the G band, around 1580 cm −1 . For samples with angles between 4.3 • and 5.9 • the ratio of the intensities of the G bands of TBG and single layer graphene ( I TBG /I SLG ) is around 2, an expected result for conditions far from resonance with the van Hove singularities (vHs) in the density of states of the TBG. On the other hand, the G band is strongly enhanced for samples with angles between 7.4 • and 10.2 • , showing in this case the occurrence of resonance with transitions between vHs in the valence and conduction bands of the TBG. The different resonance regimes shown in Fig. 1a are explained considering that the Moiré peaks in the TBG Raman spectrum can be activated by both the intralayer and interlayer resonance processes 28 , whose schematic diagrams are shown in Fig. 2a. The black and grey curves in Fig. 2a represent the Dirac cones and electronic states of the two twisted graphene layers. In the intralayer case, one incident photon with energy E intra (represented by the blue arrow) creates one electron-hole pair in one layer (in the black curve) and, then, the excited electron is scattered by a phonon with wave vector q M (represented by the red arrow) to another state of the same layer (also in the black curve). The interlayer process is represented in the right side of Fig. 2a. Now, one incident photon with energy E inter (green arrow) creates one electron-hole pair in one layer (in the black curve) and, then, the excited electron is scattered by a phonon with wave vector q M to a state of the other layer (in the grey curve). Notice that the interlayer resonance energy E inter coincides with the energy separation between the vHs ( E vHs ), depicted also as a green arrow in the left side of Fig. 2a, which describes the interlayer process for phonons with q = 0 momentum.
The dependence of the intralayer and interlayer resonance energies ( E intra and E inter ) on the twisting angle θ are shown as dashed blue and dashed green lines in Fig. 2b, respectively 28 . Therefore, using the dependence of E intra and E inter on θ and considering the excitation energies of each spectrum in Fig. 1a, we conclude that TBGs with angles in the range 4 • -6 • (highlighted by the light blue strip in Fig. 1b) are in resonance with the intralayer process, whereas TBGs with angles in the range 7.5 • -10.2 • (highlighted by the light green strip in Fig. 1b) are in resonance with the interlayer process. These conclusions are nicely in agreement with the behavior of the G band in Fig. 1a, which does not exhibit an enhancement for samples in the 4.3 • -5.9 • range, a typical behavior of the intralayer resonance process, but is hugely enhanced for samples in the 7.4 • -10.2 • range, as expected for the interlayer resonance. Figure 2b also shows that the Moiré peaks for the magic angle TBG ( θ = 1.1 • ) are expected to be enhanced in the Raman spectrum for excitation energies around 0.40 eV and 0.24 eV for the intralayer and interlayer processes, respectively, which are much smaller that the excitation energies used in this work.
Let us now discuss how the intralayer and interlayer processes enhance the distinct phonons of graphene. We can observe in Fig. 1a,b that only phonons of the TA and the LO branches appear in the case of the intralayer resonance, whereas phonons from all branches (ZA, TA, LA, oTO, iTO and LO) are enhanced by the interlayer resonance process. Figure 2c shows the Raman spectra of four different samples of TBGs with angles 4.3 • , 5.9 • , 9.2 • and 10.2 • in the region of the iTO and LO peaks recorded with the excitation energies of 1.50 eV, 1.82 eV, 2.41 eV and 2.71 eV, respectively. The excitation energy for each sample corresponds to the situation of maximum resonance observed for the Moiré peaks. Notice that the ratio of the intensities of the G band in TBG and single layer graphene ( I TBG /I SLG ) in Fig. 2c is always around two, showing thus that the Moiré peaks are enhanced by the intralayer el-ph process. Interestingly, the intensity of the LO peak at resonance condition strongly depends on the twisting angle θ . For example, in the case of the 4.3 • sample in Fig. 2c, it becomes four times more intense than the G band. For the θ = 5.9 • sample, the LO peak is not so strong but still more intense than the G band. For the 9.2 • sample the LO peak is much weaker that the G band and, in the case of the 10.2 • sample, we can only observe a weak feature at the position of the LO peak in Fig. 2c. Concerning the activation of the iTO phonons by the intralayer process, we do not observe any peak associated with the iTO phonons in the spectra of TBG  Figure 2d shows the Raman spectra of four different samples of TBGs with angles 7.1 • , 10.2 • , 14.6 • and 16.3 • in the region of the iTO and LO peaks recorded, respectively, with the excitation energies of 1.39 eV, 1.82 eV, 2.41 eV, and 2.66 eV. The excitation energy for each sample corresponds to the situation of maximum resonance observed for the Moiré peaks. We can now observe a huge enhancement of the G band of the TBGs, which becomes at resonance condition up to 80 times more intense than the G band of single layer graphene, as shown by the dashed curves in Fig. 2d. In the case of the 7.1 • and 10.2 • samples, we can clearly observe in Fig. 2d the iTO and LO peaks below and above the G band, respectively, showing that the interlayer el-ph process enhances both the iTO and LO phonons of graphene. However, the enhancement of the LO peak is weaker in the interlayer case when compared with the intralayer el-ph process shown in Fig. 2c. Notice that the LO peak cannot be observed in the spectra of the 14.6 • and 16.3 • samples, as shown in Fig. 2d. The LO frequency is closer to the G band for these angles (see Fig. 1b) and the appearance of the LO peak in the spectra may be masked by the huge enhancement of the G band.
The double-resonance Raman model for the intralayer and interlayer el-ph processes. In the last sections, we have shown that the intralayer and interlayer processes enhance phonons from different branches in a distinct way and that some phonon branches are not enhanced by the intralayer resonance process. Now we will address these points considering the double-resonance (DR) Raman model for graphene and the specific momentum conservation conditions for the intralayer and interlayer electron-phonon processes. In graphene, phonons with q = 0 can be accessed through the DR Raman process either by a defect induced mechanism or by a scattering of two phonons with opposite momenta q and −q 7,16 . The intralayer process in TBG can be considered as a DR mechanism that occurs in one graphene layer where now the momentum conservation is guaranteed by the Moiré vector q M (θ) instead of being provided by a defect or an additional phonon. We can make a parallel between the LO intralayer process in TBG with the defect-induced D ′ band in graphene, around 1620 cm −116 . The D ′ band originates from the intravalley DR Raman process in graphene involving one phonon www.nature.com/scientificreports/ of the LO branch, where momentum conservation for the process is provided by a defect. For a TBG, momentum conservation is provided by the periodic potential generated in one layer by the other one twisted by θ.
It is worth to mention that the iTO phonons are not activated by the intravalley double resonance process in disordered graphene, and this result was ascribed to the destructive quantum interference effect that suppress their contribution in the spectrum 27 . Additionally, it is known that the wave vectors in the BZ of graphene along the Ŵ − K − M direction belong to the C 2v point group and their symmetries are labeled by the letter T. The electronic bands have either T 2 or T 4 symmetries 31,32 . Two excited states in the same valley have necessarily different symmetries. Pair of excited states with same symmetry ( T 2 or T 4 ) occurs only if electrons are in different valleys, around the K and K ′ points. Along this high-symmetry direction, the iTO and LA phonons have T 1 symmetry and the LO and TA phonons have T 3 symmetry 31,32 . Group theory shows that electronic transitions between excited states of same symmetry require a T 1 -symmetry phonon, whereas a T 3 -symmetry phonon only allows transitions between excited states of different symmetries. Therefore, the intravalley DR process along the Ŵ − K − M direction, where the excited states have different symmetries, involves only LO and TA phonons, that have T 3 symmetry. The experimental results in this work for the intralayer resonance enhancement in small angle TBGs are in agreement with the group-theory predictions along the high symmetry Ŵ − K − M direction. However, the lack of observation of the iTO peaks even in small twisting angles regime cannot be explained solely by group-theory, since the direction of the q M mapped on the BZ of one graphene layer always deviate from the high-symmetry Ŵ − K line for θ = 0 (dotted white lines in panels of the Fig. 4).
In the DR model, the Raman intensity related to a phonon belonging to branch ν in the phonon dispersion is calculated by means of fourth-order time-dependent perturbation theory. Here we consider a sum over all electron wave vectors ( k ) in the SLG BZ, but we take in account only phonons with wave vectors equal to q M , for a given twisting angle.
The expression above corresponds to the Fermi golden rule transition probability for the 4-th order process 27,33 , where the system in an initial state i goes through three intermediate states A, B and C, to a final state f. The terms in the numerator corresponds to the matrix elements of the four processes: from right to left, M iA corresponds to the matrix element of the absorption of light and creation of an electron-hole pair, M AB corresponds to the matrix element of the electron-phonon (el-ph) interaction that scatters the carriers from A to B, M BC is related elastic scattering provided by the Moiré potential, and M Cf corresponds to the matrix element of the recombination of the electron-hole pair and emission of the scattered photon. ε L is the excitation energy, ε A , ε B and ε C are the energies of states A, B and C, respectively, and γ is the broadening energy related to the finite lifetime of the intermediate states. In this work M iA , M Cf , ε A , ε B , ε C and γ are calculated exactly as in Ref. 27 and M BC is considered as a constant.
The el-ph matrix element is given by M AB = �k ′ |�H q,ν |k� and it connects the electronic states |k� and |k ′ � . �H q,ν is the derivative of the Hamiltonian with respect to the periodic displacement related to the phonon from branch ν with wave vector q 27 . In the intralayer process, the electronic state |k� in one layer differs from the state |k ′ � in the same layer by the Moiré vector q M , as depicted in Fig. 2a, while in the interlayer process the difference is given by a rotation of θ , since |k� and |k ′ � are in different layers, followed by a scattering with a Moiré vector q M . In other words, we consider the conditions: and for the intralayer and interlayer processes, respectively, being R(θ) the rotation matrix. Figure 3a shows the total calculated Raman intensity obtained by the sum of intralayer and interlayer processes described above, for θ between 4.3 • and 10.2 • . We show here only the results for excitation energy of 1.50 eV, but the main results are qualitatively similar for excitation energies in the NIR and visible range. The overall agreement between calculations of Fig. 3a and the experimental results shown in Fig. 1a is clear. In particular, we can also distinguish in the calculations the two regimes where the intralayer or interlayer el-ph mechanisms are more relevant. For angles between 4.3 • and 5.3 • we see that the peaks associated with the TA and LO modes are more intense than the LA and iTO ones, respectively, in agreement with the experimental results. On the other hand, for angles between 5.9 • and 10.2 • , the TA, LA, iTO and LO modes have all non negligible intensities. All spectra were normalized by the respective LO intensity, but the acoustic (TA and LA) peaks were multiplied by a factor of 100 for a better visualisation. Since we do not consider out-of-plane polarized light, the ZA and oTO modes cannot be described in this model. We remark that the G band is a single-resonant process, which is not included in the present calculations (the vertical grey strip in Fig. 3a represents the region where the G band is supposed to appear in the spectra).
In Fig. 3b we plot separately the calculated spectra considering the contribution of the interlayer (green curves) and intralayer (blue curves) processes for TBGs with twisting angles of 4.5 • (bottom panel) and 8.5 • (upper panel), and the experimental spectra (red curves) for the same angles shown in Fig. 1a. Notice that the enhancement of the peaks and the relative intensities between the longitudinal and transverse components of the acoustic and optical phonons follow the same trend observed in the experimental Raman spectra. For the 4.5 • sample, we only observe the enhancement of the TA and LO peaks by the intralayer process whereas, for www.nature.com/scientificreports/ the 8.5 • sample, the enhancement is provided by the interlayer el-ph process and phonons of the TA, LA, iTO and LO branches are activated. These results show that, indeed, the intralayer process is stronger for small angle TBGs, while the enhancement for TBGs with intermediate angles is dominated by the interlayer el-ph process.
In order to achieve a deeper understanding regarding the contribution of each phonon branch for the calculations shown in Fig. 3a,b, we have also performed calculations using Eq. (1) considering all phonon wave vectors q , without restricting the calculation to the Moiré wave vectors ( q M ). The intensities of the DR process Ĩ ν (q) as a function of the phonon wave vector q , where ν represents the iTO, LO, TA and LA phonon branches, are shown in Fig. 4a-d for all phonon wave vectors q within the graphene BZ represented by the black hexagons. The white dotted lines in these figures represent the Moiré vector q M paths as function of the twisting angle, being θ = 0 at the Ŵ point and θ = 30 • at the middle point between K and M. We see in Fig. 4a,b that the iTO peak intensity for q-vectors close to the Ŵ point is almost null, while the LO peak is very intense in this region. For the acoustic peaks shown in Fig. 4c,d, we see that the TA and LA peaks are more intense close to Ŵ , indicating that both peaks could be seen in the Raman spectra for small θ TBGs. Furthermore, we have computed the differences between Ĩ ν (q) of longitudinal and transverse branches.  Figure 4e,f show that TA and LO are the most intense peaks in the spectra when θ is small (DR processes close to Ŵ ) but, for intermediate twisting angles, both longitudinal and transverse components of the acoustic and optical peaks have comparable intensities. We can thus conclude that the experimental results can be very well explained and understood by the double resonance (DR) Raman calculations for graphene considering momentum conservation provided by the Moiré potential and the intralayer and interlayer el-ph resonance conditions.

Conclusion
The Raman spectra of twisted bilayer graphene exhibit a number of extra peaks associated with phonons of the acoustic and optical branches of graphene with the wave vector of the Moiré pattern, which can be activated by both the intralayer and the interlayer electron-phonon processes. In this work, we studied how the intralayer and the interlayer el-ph processes enhance the intensity of each extra peak in the Raman spectrum of twisted bilayer graphene (TBG). CVD grown samples of TBGs with different twisting angles θ in the range 4 • -16 • were measured using several laser excitation energies in the range 1.39-2.71 eV. As predicted by the dependence of the All spectra were calculated using the double-resonant Raman model described in the text, which includes interlayer and intralayer el-ph processes, for the excitation energy of 1.50 eV. The intensities in the acoustic region ( ω < 700 cm −1 ), which shows the peaks associated to TA and LA modes, were enhanced by a factor of 100 for a better comparison with the optical region ( ω > 1400 cm −1 ), that shows the iTO and LO peaks. The grey strip between 1570 cm −1 and 1600 cm −1 denotes the region of the G band, which is a first-order process not included in the calculations. (b) Theoretical Raman spectra of TBG structures for θ = 8.5 • and 4.5 • showing separately the intralayer (blue curves) and interlayer (green curves) processes, for the excitation energy of 1.50 eV, compared with experimental measurements (red curves). In this case, the acoustic region was magnified by a factor of 50 compared to the optical one and all the theoretical frequencies shown in the figures were shifted by ±20 cm −1 , in the optical and the acoustic region, respectively. Peaks associated with all phonon branches of graphene were observed experimentally in the Raman spectra of TBG samples in resonance with the interlayer process. They can be as intense as the G band of single layer graphene, but much weaker than the TBG G band. Their positions are nicely described by the calculated dependence of phonon frequencies on θ . On the other hand, our experimental results shows that the intralayer process enhances preferentially phonons of the TA and LO branches, considering the laser energies and the range of θ of the samples used in this work. The Raman enhancement by the intralayer process is very intense and the Moiré peak can be much stronger than the G band at resonance condition. It was observed that the enhancement of the LO peak increases with decreasing twisting angle θ . The experimental results are explained by theoretical calculations of the double-resonance Raman process in graphene, within the fourth-order perturbation theory approach, considering the momentum conservation provided by the Moiré vector q M (θ) for both the intralayer and interlayer el-ph processes. The observed enhancement of the acoustic and optical phonons by the intralayer and interlayer processes was nicely reproduced in the calculated DR spectra, showing the importance of taking into account the symmetry-related selection rules for the electron-phonon interaction in graphene, as well as the quantum interference that affect the intensities of the Raman peaks activated by the Moiré superlattice.

Methods
Raman measurements of the TBG samples with different twisting angles θ were performed using several laser lines in the near infrared (NIR) and visible ranges. Experiments in the visible range were performed in a DILOR XY triple-monochromator spectrometer equipped with a N2-cooled charge-couple device detectors, using an Ar/Kr laser with 12 excitation energies in the visible range. The NIR Raman measurements were performed on a home-made setup, including an iHR-550 Horiba spectrometer equipped with a liquid-nitrogen-cooled silicon CCD detector. A tunable CW Ti:sapphire laser filtered using tunable laser line filters 34 was used for excitation. The scattered light was collected through a × 100 objective (NA = 0.95) using a backscattering configuration.
The TBG samples were grown by CVD technique using methane (99.99 per cent) on polycrystalline Cu foils 30 . The graphene layers were first covered by a thin layer of polycarbonate, followed by etching in HCl aqueous solution to remove the Cu in the transfer process. The polycarbonate film with attached graphene was then transferred onto fused silica. Finally, the polycarbonate film was removed using chloroform.
The double-resonance (DR) Raman intensities were computed in the fourth order perturbation theory. The eigenenergies and the eigenstates were obtained through a fifth neighbour tight-binding model fitted to a DFT + GW calculations, as well as the phonon dispersion used in all calculations. The electronic wave vectors www.nature.com/scientificreports/ were distributed in a 480 × 480 k-point grid over a graphene BZ and only phonon wave vectors q which conserve momentum with the Moiré vector q M (θ) for a given twisting angle θ were used in TBG, while in the full DR Raman calculation in SLG, a 240 × 240 q-point grid for all six phonon branches were used. Despite the restrictions associated to the interlayer and intralayer processes, all Raman calculations were performed following the procedures described in Ref. 27 .